/**
 * Copyright 2012 Brigham Young University
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *    http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package edu.byu.nlp.math;

import org.apache.commons.math3.special.Gamma;

import edu.byu.nlp.util.DoubleArrays;

/**
 * @author rah67
 *
 */
public class GammaFunctions {

	private static double EPS = 1e-10;
	private static double FALLING_THRESHOLD = 8;
	
	private GammaFunctions() { }
	
	public static double logRatioOfGammas(double numerator, double denominator) {
		// There are speed ups when the difference is integral
		double diff = numerator - denominator;
		if (Math2.isIntegral(diff, EPS)) {
			return logRatioOfGammasByDifference(denominator, (int)Math.round(diff));
		}
		
		return Gamma.logGamma(numerator) - Gamma.logGamma(denominator);
	}
	
	public static double logFallingFactorial(double x, int k) {
		// A few optimized special cases
		if (k == 0) {
			return 0.0;
		}
		
		if (k == 1) {
			return Math.log(x);
		}
		
		double acc = 0.0;
		for (int i = 0; i < k; i++) {
			acc += Math.log(x + i);
		}
		return acc;
	}

	public static double logRatioOfGammasByDifference(double x, int diff) {
		// If diff is negative, then invert the fraction (in log space)
		if (diff < 0) {
			return -logRatioOfGammasByDifferenceUnchecked(x, -diff);
		}
		
		return logRatioOfGammasByDifferenceUnchecked(x, diff);
	}

	private static double logRatioOfGammasByDifferenceUnchecked(double x, int diff) {
		assert(diff > 0);
		// We can use a falling factorial when the difference is small.
		if (diff < FALLING_THRESHOLD) {
			return logFallingFactorial(x, diff);
		}
		
		return Gamma.logGamma(x + diff) - Gamma.logGamma(x);
	}

	public static double logBetaSymmetric(double alpha, int dim) {
		return dim * Gamma.logGamma(alpha) - Gamma.logGamma(dim * alpha);
	}

	public static double logBeta(double[] alpha) {
		double sum = 0.0;
		for (double a : alpha) {
			sum += Gamma.logGamma(a);
		}
		return sum - Gamma.logGamma(DoubleArrays.sum(alpha));
	}
}
